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Abstract 

This article is devoted to the numerical study of various finite difference approximations 
to the stochastic Burgers equation. Of particular interest in the one-dimensional case is 
the situation where the driving noise is white both in space and in time. We demonstrate 
that in this case, different finite difference schemes converge to different limiting processes 
as the mesh size tends to zero. A theoretical explanation of this phenomenon is given and 
we formulate a number of conjectures for more general classes of equations, supported by 
numerical evidence. 
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1 Introduction 

This article studies several finite difference schemes for the viscous stochastic Burgers equation: 
d t u(x,t) = vd*u-g(u)d x u + a£(x,t), a;e[0,27r], t > 0. (1) 



In this equation, £ denotes space-time white noise, that is the centred, distribution-valued Gaus- 
sian random variable such that E(£(x, t)£(y, s)) = S(t — s)5(x — y). We will always endow this 
equation with periodic boundary conditions and we will consider solutions u taking values cither 
^- ■ in R or in M" (in which case g is matrix- valued in general). 

Motivations for studying the stochastic Burgers equation are manifold. Just to name a few, 
it is used to model vortex lines in high-temperature superconductors [BFG + 94j . dislocations in 
disordered solids and kinetic roughening of interfaces in epitaxial growth [Bar96 , formation of 
large-scale structures in the universe [GSS85, SZ89 , constructive quantum field theory [BCJL94 , 
etc. Since in the case a = and g(u) oc u this equation is furthermore explicitly solvable via 



in 

the Hopf-Cole transform (u = d x v/v, where v solves the heat equation), it comes as no surprise 
that a wealth of numerical and analytical results are available. From a purely mathematical 
point of view, let us mention for example the well-posedness results from [BCF91. BC JL94J 
IDPDT941 |Gyo98[ lKim06] and the ergodicity results obtained in }TZ06[ IGM05] . One remarkable 
achievement was the construction of a stationary solution in the inviscid limit with non- vanishing 
noise EKMS00 (dissipation then occurs purely through shocks). From a more quantitative 
perspective, the scaling exponents of the solutions in the small viscosity limit have attracted 
considerable interest, both in the physics and the applied mathematics literature [YC96, BM96, 
IKra99l lEVEOOal lEVEOOb] . 

The white noise term £ in (JT|) leads to solutions u which, in general, will be very "rough" . 
In particular, u as a function of the spatial variable x will not be differentiable in the classical 
sense, but will only possess some Holder regularity. As a consequence, it transpires that the 
solutions to ([1]) are extremely unstable under natural approximations of the nonlinearity, and 
this is the phenomenon that will be explored in this article. For example, for any a, b > with 
a + b > 0, one can consider the approximating equation 

d t u £ (x, t) = v d 2 x u £ - g(u e ) D £ u £ (x, t) + a £(x, t), (2) 

where the approximate derivative D e is defined as 

u{x + ae,t)-u{x-be,t) 
D £ u(x,t) = . 3 

(a + b)e 



In the absence of the noise term £ it would be a standard exercise in numerical analysis to show 
that the solution of ^ converges to the solution of ([T]) as e — > 0. This is just an example of the 
widely accepted 'folklore' fact that.if an equation is well-posed, any 'reasonable' approximation 
will converge to the exact solutionis 

In this article, we will argue that, if £ is taken to be space-time white noise, the limit of 
@ as £ — > depends on the values a and b and is equal to ([T]) only if either g is constant or 
a = b\ Furthermore, it will follow from the argument that, if one considers driving noise that is 
slightly rougher than space-time white noise (taking a noise term equal to (1 — d x ) a dw(t) with 
a £ (0, 1/4) still yields a well-posed equation), one does not expect solutions to the approximate 
equation ([2]) to converge to anything at all, unless a = b. Our methodology here is to first 
present a heuristic argument which allows to derive quantitative predictions for the effect of the 
finite difference discretisation on the solution. We will then use numerical experiments to verify 
these predictions. 

At this point we would like to emphasise that the aim of this article is certainly not to 
advocate the use of a finite difference scheme of the type <j2j) to effectively simulate (JXJ) - Indeed, 
we will show in section T2.2I below that approximations of the nonlinearity of the type D e G(u), 
where G is the antiderivative of g already have much better stability properties. Instead, our 
aim is merely to give a striking illustration of the fact that caution should be exercised in the 
simulation of stochastic PDEs driven by spatially rough noise. 

The text is structured as follows: We start in section [2] by presenting our argument for the 
case of the stochastic Burgers equation, i.e. for g{u) = u. In section [3] we will present the 
corresponding results for more general equations and in section 2] we study the limit of vanishing 
noise and viscosity. Finally, in appendix [A] we discuss some technical aspects of the simulations 
used throughout this article. 
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2 Stochastic Burgers Equation 

In this section we consider the stochastic Burgers equation 

du = v d x u dt — u d x u dt + a dw, (4) 

as well as the approximation given by 

du e = v d 2 x u E dt - u e D e u £ dt + a dw. (5) 

The solutions u and u e take values in the space £ 2 ([0, 2tt],K) on which the operator d x is 
endowed with periodic boundary conditions, and w is an L 2 -cylindrical Wiener process (see 
[DPZ92] for details). Equation (£[]) is well-posed, since we can rewrite ud x u as ^d x (u 2 ) which 
is locally Lipschitz from the Sobolev space H 1 / 4 into H' 1 , thus allowing us to apply general 
local well-posedness theorems as in [DPZ92, Hai09 . For fixed time t > 0, the solutions to (QJ 
have the regularity of Brownian motion when viewed as a function of the spatial variable x. In 
particular, they are not differentiable in x. Figure [T] shows numerical solutions of equation (O 
for different values of a and b. One can see that different choices for these parameters lead to 
an 0(1) difference in the solutions. 

Our aim in this article is to understand and quantify these differences. In particular, in 
conjecture [1] below, we compute a correction term to (UJ) and we verify numerically that the 
solutions to ([5]) converge to the corrected equation as e — > 0. This understanding will then allow 
us to conjecture the appearance of similar correction terms in more complicated situations and 
we will again verify these conjectures numerically. 

1 Remember that we are working at fixed non-zero viscosity here, so there is no ambiguity in the concept of 
solution and we do not require an upwind scheme in order to obtain convergence. 




Figure 1: Illustration of the discretisation error for the finite difference method f$). The three 
lines correspond to a right-sided discretisation (a — 1, b — 0; top-most curve), a centred discreti- 
sation (a = 1. b = 1; middle curve) and a left-sided discretisation (a = 0. b = 1; bottom-most 
curve), all computed using the same instance of the driving noise. The picture clearly shows that 
there is an 0(1) difference between the solutions obtained by the three different discretisation 
schemes. From the argument presented in the text, we assume that the exact solution of Op will 
be closest to the middle of the three lines. 



2.1 Heuristic Explanation 

For simplicity, instead of u we consider the solution v to the stochastic heat equation 

dv = v d x v dt + a dw(t). (6) 

Since the properties of the discretisation of differential operators only depend on local properties, 
and since v has the same spatial regularity as u, it will be sufficient to study how well v D e v 
approximates vd x v = ^d x v 2 . 

By expressing v in the Fourier basis {e m:E /v / 27r} _- it is easy to check that the stationary 
solution to © is 



v(t,x) 
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where the £o is a (real-valued) standard Brownian motion and £„ for ti^O are complex-valued 
Ornstein-Uhlenbeck processes with variance 1 (in the sense that E|£„(£)| 2 = 1) and time constant 
un 2 which are independent, except for the condition that £_ n = £„. Therefore, the derivative of 
v is given (at least formally) by 



g£n(f)e* 



d x v(x) = ^ 
The e-approximation to the derivative given in (0) , on the other hand, is given by 



D e v{x) = J2 



a£ n (t)e 
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(7) 



(8) 



It is clear that the terms in ((5J) are a good approximation to the terms in (|7|) only up to n ~ e . 
For larger n, the multiplier in ((5J) will decrease like n . 

For our analysis we restrict ourselves to the constant (n = 0) Fourier mode. Our numerical 
experiments, below, show that the contributions from this mode are already enough to explain 
the observed differences between the solutions of the approximating equation (O and the exact 



solution. Since v d x v is a total derivative, the 0-mode of this term vanishes. In contrast, the 
0-mode of v D F v does not vanish at all: We obtain instead for this term the sum 



1 . 1 v^ cr 2 4_ n (i)C„(t) e* 
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and the expectation of this expression, as e — > 0, converges to 
\imE(-vD 6 v)(x) = \imE(-vD 6 v | -7=)- 
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As a consequence, one expects the following result. 

Conjecture 1 The solution of the approximating equation {3J) converges, as e — > 0, to t/ie 
solution of 

du = v d^u dt — u d x u dt H dt + a dw. (10) 

4i/ a + 

Thus, the approximation converges to the stochastic Burgers equation Ol) only for a = b. 

Remark 2.1 The solution to the stochastic Burgers equation (or rather the integrated process 
which solves the corresponding KPZ equation) arising as the fluctuations process in the weakly 
asymmetric exclusion process BG97, BQS09 is driven by the derivative of space-time white 
noise. As a consequence, it does not solve an SPDE that is well-posed in the classical sense and 
can currently only be defined via the Hopf-Cole transform. Such a process is even rougher (by 
"one derivative") than the process considered here and one would expect the 'wrong' numerical 
approximation schemes to fail there in an even more spectacular way. 

Remark 2.2 One may think of two reasons why the correction term vanishes when a = b. On 
the one hand, this discretisation is more symmetric. On the other hand, it yields a second-order 
approximation to the derivative at x. The correct explanation is closer to the first one. Indeed, 
consider for example the discretisation 

t^ \, \ cu(x + 2e) + (1 - 3c) u{x + e) + 3cu(x) - (1 + c) u(x — e) . 

(D s u){x)tt — . (11) 

This discretisation is of second order for every c e R. However, if we perform the same calculation 
as above with this discretisation, we obtain a correction term equal to 

a 2 f°° c cos 2k — 4c cos k + 3c ca 2 

2tw J q 2k 2 K_ ~~8i7' 

which vanishes only if c = 0, i.e. when dill) conincides with the symmetric discretisation (|3"]l. 

In the above calculation, both the limiting equation and the approximating equation live in 
the same space. It is possible to carry out a similar analysis in the case where the approximating 
equation takes values in a different space, typically R N for some large N. For example, we can 
consider the 'finite difference' approximation given by 

„, u(x + 5) - 2u(x) + u(x — S) d «f / A \ , s 
d x u*i = (A N u)(x) 

(12) 
„ n , <u(x + 5)-u(x) dcf 

uo x u w uDgu — u(x) j = i<n{u){x), 



where we set 5 = -£■ for the approximation with N gridpoints. In this setting, we approximate 
u by u N € R^ with u 1 ^ ~ u(JS). The natural candidate for the approximation of space-time 



white noise is then given by dWf* = <5 -1 / 2 dWj, where the Wj are independent, one-dimensional 
standard Brownian motions. This is the correct scaling, since it ensures for example that for every 
continuous function g, ^ W^g(5j) is a Wiener process with covariance J^. g(Sj) 2 5 s=s J g 2 dx. 
With this notation, we consider the approximating equation given by 

du N = vA N u N dt - F N {u N ) dt + <rdW N {t). (13) 



Let us take N even for the sake of simplicity. It is then straightforward to check that the 

N . i JV 

2 T J- 1 • • ■ > 2 



eigenvectors for A at are given as before by e mx with n. = —^ + 1,...,^, but the corresponding 



eigenvalues are 

, 2 , . /2 . /W(5s\ 2 dof 2 

A ™ = ^2( cosn(5 " 1 ) = ~(j slI Hy)J = ~ ?7 ™- 

(Note that for fixed n and small 6, one has indeed A„ w — ?i 2 .) It then follows as previously that 
the solution to the linearised equation is given by 

and that its discrete derivative -D^v is given by 

-^ ae inx Z„{t) c inS - 1 



D 5 v(t,x) 






Note that the sums in both expressions only run over the values n = — — + 1 , . . . , -j . Similarly 
to above, we obtain that the expectation of the zero mode of the product v D$v is given by 

1 1 2^/2^1 2 N / 2 r 2 

, | 1 > 1 cr z v ^cosdn-l (T z v-^ d cr z 

v v n=l '" n=\ 

One therefore expects the following result. 

Conjecture 2 27ie solution u N of the finite difference approximation US]) converges, as N — > 
oo, to the solution of 

2 

du = v d-rUdt — u d x u dt H dt + a dw(t) . 

Au 

To test this conjecture, we use the following numerical experiment: We numerically solve 
both the "approximating" equation (fT3"|) and the "corrected" SPDE 



(T 
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du 7 = vd x u-y dt — Uj d x u-y dt + 7 — dt + a dw(t), (15) 



until a fixed time T, using the same instance of the noise. To discretise the term ud x u in (fH 
we use the approximation (u(x) 2 /2 — u(x — S) 2 /2J/5, which is known to converge to the exact 
solution, see section 12.21 below. Also, for increased accuracy, we use a finer grid for (fT5|) than 
we did for p3[) . Now we can compare the two solutions by considering || u (T, ■ ) — u 7 (T, - ) |] 2 
as a function of 7: If conjecture [5] is correct, this function will take its minimum at 7 ~ 1/4. 
The result of such a simulation is given by the line labelled "finite difference" in figure O The 
minimum of the curve is indeed located close to 7 = 1/4, thus giving support to conjecture [5J 

Even though the correction terms in conjectures [T] and [2] (with a — 1 and b — 0) coincide, the 
constants arise in completely different ways. This might lead to the speculation that the value of 
this constant is an intrinsic property of the limiting equation, rather than of the particular way 
of approximating it. The following argument shows that this is not the case. One can repeat the 
calculation leading to conjecture [5] with a 'spectral Galerkin' approximation of the linear part of 
the equation, but retaining a 'finite difference' approximation of the nonlinearity. In other words, 
we consider (fl3|) as before, but we take for Ajv the self-adjoint matrix with eigenvectors e mx and 
eigenvalues —n 2 . (This can be achieved by first applying the discrete Fourier transform, then 
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Figure 2: This figure compares the solution u N of the finite difference approximation US]) to the 
solution u 7 of the "corrected" SPDE \15\) (which includes an additional drift term ja 2 jv). The 
curves show the L 2 -norm difference between the solutions (for the same instance of the noise) as 
a function of 7, once using finite difference discretisation ilty) of the linear part (corresponding 
to the top-most curve in figure [7J) and once using the spectral Galerkin discretisation. The two 
vertical line segments give the predicted locations for the minima of the two curves. It can be 
seen that predictions and simulations are in good agreement. 



multiplying the nth component by —n 2 , and finally applying the inverse Fourier transform.) In 
this case one has r\ n = n and it transpires that the correction term is given by 
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which is clearly different from ([Fiji . 

To verify that the spectral Galerkin discretisation of the linear part indeed leads to this 
different correction term, we modify the code which we used to validate conjecture [5] above. The 
result of this simulation is given by the line labelled "Galerkin" in figure [2 Again, there is good 
agreement between our conjecture and the simulation results. 



2.2 The Case of More Regular Noise 

To conclude this section, let us argue that the situation considered in this article is truly a 
borderline case in terms of regularity and that if we drive ([5]) by noise that gives rise to slightly 
more regular solutions, the solutions of (J5J) converge to those of ([U without any correction term. 
Indeed, consider a general semilinear stochastic PDE driven by additive noise: 



du = -Audi + F(u) dt + Q dw(t), 



(16) 



where A is a strictly positive-definite selfadjoint operator on some Hilbert space H, F is a 
(possibly unbounded) nonlinear operator from H to T-L, IF is a standard cylindrical Wiener 
process on H, and Q : % — > % is a bounded operator. 

Denote furthermore H a = V(A a ) for a > and let %~ a be the dual space to W a under the 
dual pairing given by the Hilbert space structure of H. (So that H~ a is a superspace of H for 
a > 0.) We also denote by || • || a the natural norm of "H a . Finally, we denote as before by v the 
solution to the linearised system 

dv = —Av dt + Q dw(t), 

which we assume to be an "H-valued Gaussian process with almost surely continuous sample 
paths. One then has the following convergence result: 



Theorem 2.3 Assume that there exists 7 > and a £ [0, 1) such that F: IF 1 — > W~ ( ~ a is locally 
Lipschitz continuous and such that the process v has continuous sample paths with values in TF 1 . 
Assume furthermore that F £ : T-L 1 — > r hF~ a is such that F s is locally Lipschitz and such that the 
convergence F e (u) — > F(u) takes place in 'H' y ~~ a , locally uniformly in TF 1 . Then, the solutions u e 
to 

du e = -Au e dt + F(u £ )dt + Qdw(t), (17) 

converge to the solutions to H16\) as e — > 0. 

Proof. The proof is straightforward and we only sketch it. We assume without loss of generality 

that the parameter e is chosen in such a way that for every R > there exists a constant Cr 

such that 

sup \\F e (u) - F(u)|| 7 _ a < C R B. (18) 

IMI-r<-R 

Denote now by u the solution to (I16[) with initial condition uq and by u e the solution to (1171) with 
initial condition Uq. Let t > be small enough so that ||u(s)|| 7 < R and \\u(s) — u e (s)\\ 1 < R 
for s < t. We then have 

\\u(t) - u s (t)h < ll«o - uoll-y + C f (t- s)- a \\F £ (u e (s)) - F(u(s))\\ 7 - a ds 

Jo 

< ||«o - «oll 7 + C [ (* - s y a \\ F s(Ms)) - F(u e (*))|| 7 _ a ds 
Jo 

+ C f (t - s)- a \\F(u e (s)) - F(u(s))\\^ a ds 



11 



< ||u - W0II7 + C R e + Cut 1 a sup \\u(s) - « e («)|| 



t-1 — a , 

l7~ 

8<t 



The claim then follows by taking t sufficiently small and performing a simple iteration. □ 

Despite its simplicity, this criterion is surprisingly sharp. Indeed, we argue that if we consider 
(U)) but with the space-time white noise dw replaced by (1 — d^)~ d dw for 5 > 0, then the 
assumptions of theorem l2.3l can be satisfied with some choice of exponent 7 for the approximation 



F e (u)(a:) - u(x) u( - X + 6) U{x) =u(x)(D e u)(x). 

Obviously, this cannot be the case when 6 = 0, since we then observe the convergence to solutions 
to mjj). 

Indeed, we first note that since the linear operator appearing in (J3J) is of second order, we 
have the correspondence 

W =H 2 '>, 

between interpolation spaces and fractional Sobolev spaces. In order to keep our notation coher- 
ent throughout this section, we still denote by || • || 7 the norm of TF 1 , i.e. ||u|| 7 = ||(1 — c^) 7 m||£2, 
where we implicitly endow d% with periodic boundary conditions. With this notation, we get 
the following result. 

Lemma 2.4 For 7 > and a £ [-|, 1], we have \\D E u — 9j;u|| 7 _ a < Ce 2a ~ 1 ||zi|| 7 . 

Proof. The operator D e — d x is given by the Fourier multiplier M £ (k) = e _1 (e lke — 1 — ike). We 
immediately obtain the bound 

|M 8 (fc)| < k(ekM) < k{ek) 2a ~ l , 

from which the claim follows at once. □ 

Since furthermore TF is an algebra for 7 > j, we conclude that the bound (fT8|) holds (for 
some different power of s), provided that we choose 7 £ (j, |] and a = 27. On the other hand, 
the solution to the linearised equation 

dv = d 2 x v dt + (1 - d 2 x )- s dw 



<J = -l/4 8=-l, 



40 
35 
30 ft 
25 | 
20 f 
15 [- 
10 h 

5 





0.0 



0.2 



0.2 



<y=— i/i6 



-\ 1 — 

= 1,6 = 1 

a=l, b=0 



0.4 



0.2 



j_ 



0.4 
t 



0.6 



0.8 



Figure 3: Illustration of the divergence of the right-sided discretisation for noise rougher than 
space-time white noise. The three panels show, for different values of 8, the L 2 -norm of numerical 
solutions to equation i20\) . The full lines correspond to a centred discretisation (a — 1, b = 1) 
whereas the dotted lines correspond to the right- sided discretisation (a — 1, b = 0). The figures 
show that only the centred discretisation seems to be stable for 8 < 0. 



belongs to 7-L 1 if and only if 7 < j+ 8, thus supporting our claim that the case 8 
borderline for the applicability of theorem 12.31 

If on the other hand we make the more natural choice 

F e (u) = ±D e (u 2 ), 



is precisely 



(19) 



then it turns out that we can apply theorem 12.31 even in the case 8 — 0. Indeed, it follows 
from standard Sobolev theory (see for example |Hai09j ) that if 7 <E (§, t), the map u h-> u 2 



is locally Lipschitz from 1C into V. 13 provided that /3 < 27 — j. As a consequence, it follows 
from lemma [23 that the approximation F e given by (|T9"f converges to ud x u in the sense of (fT5)) . 



provided that 7 € (i, \) and a e 



7,1)- 



Remark 2.5 Some ad hoc numerical scheme was shown to converge to the exact solution in 
AG06_. Some other schemes are shown to converge in |GKN02j . but only in the case 8 > 
of more regular noise. A particle approximation to a specific modification was constructed in 
[GD04J . 



2.3 Numerical Verification of the Borderline Case 

We have performed numerical simulations that corroborate the argument presented in the pre- 
vious section and show that 8 = truly is the borderline case for the limiting equation to be 
independent on the type of discretisation performed on the nonlinearity. In the case 8 < 
(i.e. the case where the driving noise is rougher than space-time white noise), our preceding 
discussion suggests that the centred discretisation should converge to the correct solution for \8\ 
sufficiently small, but that the solutions to both the left-sided and the right-sided discretisations 
should diverge as e — > 0. 

In order to verify this effect, we numerically solve the SPDE 



du = v d x u dt — u d x u dt + (1 — d x 



' dw, 



(20) 



using different discretisation schemes and different values for 8. The results are shown in figure[3l 
The simulations are in good agreement with the argument outlined above. 



3 Possible Generalisations of the Argument 

In this section, we discuss a number of possible extensions of these results to more general 
Burgers-type equations. We restrict ourselves in this discussion to the case a = 1, b = 0, i.e. 
to right-sided discretisations. This is purely for notational convenience, and one would expect 
similar correction terms to appear for arbitrary values of a and b, just as before. 

3.1 More General Nonlinearities 

Consider the equation 

d 
dui = v d x Ui dt +y djhi(u)d x Uj dt + a duii, (21) 

i=i 

for an Revalued process u and a smooth function h : R d — ► R d with bounded second and third 
derivatives. Rewriting the nonlinearity as d x [hi(u)), we see that this equation is globally well- 
posed. As before, we consider the approximating equation 

d 
du e t (x,t) = vd 2 x ul(x,t)dt + ^2d ] h l (u £ (x,t))D e u 6 3 (x,t)dt + adw t {t). (22) 

The idea now is to introduce a cut-off frequency N and to write u 6 — u e + u e , where u e is 
the projection of u £ onto Fourier modes with \k\ < N. Since the linear part of the equation 
dominates the nonlinearity at high frequencies, one expects u e to be well approximated by v, 
the projection of v onto the high frequencies. This on the other hand is small in the L°° norm 
(it decreases like N~ s for every s < |), so that 

d 

fe=i 

It now follows from the same argument as before that the term d%hi(u £ )ikD e Vj is expected to 
yield a non- vanishing contribution for k = j in the limit e —¥ and AT — > oo. Provided that we 
keep iV <C -, this contribution will again be described by (|9]), so that we expect the following 
behaviour: 

Conjecture 3 The solution of 123)[) converges, as e — > 0, to the solution of 

2 

dui — v d x Ui dt + y^ I djhi(u)d x Uj — —djjhi(u) I dt + a dun. (23) 

3 

In the one-dimensional case we can recover conjecture [1] (for a — 1, b — 0) from conjecture [3] 
by choosing h(u) — —u 2 /2 in (|23p. 

We perform the following numerical experiment to validate the functional form of the correc- 
tion term given in conjecture [3] We numerically solve both the "approximating" equation (|22j) 
and, for p: R -^ R, the "corrected" SPDE 

2 

du p = v d 2 u p dt + h! (u p )d x Up dt — —p(u p ) dt + a dw (24) 



until a fixed time T, using the same instance of the noise. As for the discretisation of (fT5|) 
above, we use the approximation (h(x) — h(x — S))/S for the term h'(u)d x u in the proposed 
limit (|2"4")1 and we also solve (|2"4")1 on a finer grid than we use for (|2"2"j) . Finally, we numerically 
optimise the correction term p (using some parametric form) in order to minimise the distance 
Hu^T, • ) — u p (T, ■ )||2- If the conjecture is correct, we expect the minimum to be attained for 
a function p which is close to the predicted correction term h" . The result of a simulation is 
shown in figure |U The figure shows that there is indeed a good fit between the conjectured and 
numerically determined correction terms. 





0.4 


1 — i 


1 1 1 


1 




0.2 






.•••'' - 


1? 


0.0 






*-v 




-0.2 






\ 




-0.4 


i 


i i i 


i i 





-0.4 


^-^ 


0.6 


T— ( 

d 


0.5 


•< 

3 


0.4 


t+H 




o 


o.:-s 


s 






0.2 


Sn 




o 


1 


W 








-d 


0.0 



n 



-2 



-1 





u 



Figure 4: Illustration of the convergence of \22]) to 123\) for the one- dimensional example h'{x) — 
sin(x) 2 . For the figure we numerically compute the finite difference solution u N to &22fy . We 
then compute solutions u p to {21$ , using a fifth-order polynomial for the correction term p. This 
polynomial is then numerically fitted to minimise \\ u N (l, -) — u p (l, • )||2- The top panel shows the 
resulting fitted correction term —ap/4v (full line) together with the correction term —o~ 2 h"(u)/Av 
predicted in conjecture [3] (dotted line). To give an idea which range of the correction term is 
actually used in the computation, the lower panel shows the histogram of the values of u N (the 
vertical bars indicate the 5% and 95% quantiles). Between these quantiles, the graph shows a 
good fit between the numerically determined and conjectured correction terms. 
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3.2 Classically Ill-Posed Equations 

Pushing further the class of equations considered in the previous subsection, one may want to 
consider equations of the form 

d 
dui — v d x Ui dt + y^ Gij (u)d x Uj dt + f(u) dt + a dwi (25) 

for some functions G: R d — > R dxd and /: R d — > R d . If we do not assume that G has an 
antiderivative and since solutions are only expected to be a-Hbldcr continuous in space for 
a < i , it is no longer even clear what it means to be a solution to this equation. So, at least 
classically, (J25D is ill-posed and the mere concept of solutions to such an evolution equation is 
difficult to establish. 

However, the discretised equation does of course still make sense for any fixed value of s. 
Furthermore, we observed numerically that there seems to be no instability as e — ¥ 0; indeed one 
observes pathwise convergence to a limiting process. By analogy with the behaviour observed for 
the situations where (1231) is classically well-posed (i.e. when G has an antiderivative), one would 
then be tempted to define solutions to ([2"5jl to be those processes that can be obtained as limits 
as e — > of the solutions to the equation where d x Uj is replaced by its symmetric discretisation. 
In analogy with conjecture G2 we would then expect solutions of the right-sided discretisation to 
converge, as e ~ > 0, to solutions of the corrected equation 

2 

dui — i/d x Ui dt + 2^{Gij(u)d x Uj — —djGij(uj) dt + f(u) dt + a dwi. 

3 

This reasoning leads to: 

Conjecture 4 As e — > 0, the equations 

d 

dui = v d^ul dt + J2G lj (u £ )D 1 e >°u £ j dt + f(u £ )dt + adw u (26) 

3=1 

where D^' denotes the right-sided discretisation, and 

d 2 

dui ^vdlu E l dt + ^(G lJ {u e )Dl' 1 u E J - ^-djGijiu 6 )) dt + f(u 6 ) dt + a dw„ (27) 

i=i 

where D^ 1 denotes centred discretisation, converge to the same limit. 

To test conjecture 01 we use the following numerical experiment: as an example we consider 
the SPDE 

dfU =—d 2 u+—( ° cos( U2 ) - sin( Ul )\ 

o tU - ^ 2 a x u + ^ ^ gin _ Q o x u 

(28) 



cr 



4_f sm(ui)co S (ui)\ ^^ gtw 



■ cos(w2) sin(u2) 



This SPDE is of the form (|2"5")) where G has no antiderivative. SPDEs like (|2"8")l occur in the 
problem described in |HSV07[ section 9] where we argue that the stationary distribution of this 
SPDE on L 2 ([0, 27r], IR 2 ) (when equipped with appropriate boundary conditions) coincides with 
the distribution of the stochastic differential equation 

For our experiment we numerically solve the SPDEs ([2"B)) and (12"?]) and compare the solutions. 
The result is displayed in figure [51 As can be seen from the figure, the simulation results are in 
good agreement with conjecture [4j 
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Figure 5: Illustration of the convergence in conjecture^ The left-hand panel shows a numerical 
solution of the M 2 -valued SPDE A28\) at a fixed time t. Since we use periodic boundary conditions, 
the plotted graph of u(t, • ) forms a loop in M 2 . The black line in this plot was obtained using a 
centred discretisation, whereas the grey line in the background was obtained using a right-sided 
discretisation. As in figure]]] there is an 0(1) difference between the two discretisation schemes. 
The two plots on the right show the differences u\ — u\ (upper panel) and u\ — u\ (lower panel) 
between the discretisation schemes A26\) and {21^ from conjecture^ (full lines). For comparison, 
the plots also show the differences u\ — u\ and u\ — u\ where u £ is the solution of the uncorrected 
SPDE H28\) using a centred discretisation (dotted lines). The graphs support conjecture R] by 
showing good agreement between the solutions of H26]) and \21ty. 
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3.3 Multiplicative Noise 

We conclude this section by considering the equation 

du — vd 2 udt + g(u)d x udt + f(u)dw, (29) 

where g is as before and / is a smooth bounded function with bounded derivatives of all orders. 
Such an equation is well-posed if the stochastic integral is interpreted in the Ito sense |GN99j . 
(Note that it is not well-posed if the stochastic integral is interpreted in the Stratonovich sense. 
This follows from the fact that, at least formally, the Ito correction term is infinite when / is 
not constant.) 

In such a case, the local quadratic variation of the solution is expected to be proportional 
to / 2 (u), so that one expects the right-sided discretisation to exhibit a correction term propor- 
tional to g'(u)f 2 (u). More precisely, in analogy to conjecture [31 one would expect the following 
statement to hold. 

Conjecture 5 The solution of 

du = v d x u dt + g(u)D 6 u dt + f(u) dw (30) 

converges, as e — > 0, to the solution of 

du = vd 2 x udt + g(u)d x udt- —g'(u)f(u)dt + f(u)dw. (31) 

To test this conjecture we perform a numerical experiment, similar to the one for conjecture^ 
the result is shown in figure [6] The fit between predicted and numerically determined correction 
term in figure [5] is worse than in figure 0] and thus the numerical test is not entirely conclusive. 

One possible reason is that the spatial resolution of our numerical simulations may not 
be sufficient. Indeed, the argument of the previous sections is based on a spatial averaging 
of the small-scale fluctuations of the process. In the case of multiplicative noise, these small- 
scale fluctuations are themselves multiplied by a the process f(u), which is spatially quite rough. 
Therefore, this spatial averaging will hold only on extremely small scales, where f(u) is essentially 
constant. In order to be seen by the numerical simulation, these scales still need to be resolved 
at sufficient precision to have some version of the law of large numbers. 

4 Small Noise/Viscosity Limit 

One regime that is of particular interest is the small noise/small viscosity limit. If one takes 
v ex a 2 in conjectures Q] and [21 one obtains a non- vanishing correction term even for arbitrarily 
small v and o~\ It is therefore of interest to study approximations to 

du = ed 2 udt — ud x udt + y/sdw (32) 

for e < 1. 

It is well-known that, in the limiting case e = 0, finite difference schemes for the Burgers 
equations can only be used with extreme caution due to the presence of shocks in the solution. 
These shocks are jump discontinuities of the solution. The values u~,u + to the left /right of 
the jumps satisfy u~ > u + ; in other words, the jumps are always downwards jumps. For 
viscosity solutionto to the inviscid Burgers equation, shocks move through the system at velocity 
\{u + + u-). 

What happens at the formation of a shock? If the limiting non- viscous Burgers equation is 
discretised as 

o t u n = -u n ■= , (33) 

we expect the discretisation scheme to be stable only when the discretisation is upwind in the 
sense that the direction of the discretisation coincides with the direction of propagation of the 



2 Also called entropy solutions, these are the solutions that are obtained as limits of (1321) as e — > 
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Figure 6: Illustration of the convergence of i3U\) to i31\) in conjecture^ for the case g{u) = —u 
and f(u) = 1 + 5 cos(3u). See figureU\for an explanation of the graphs. Here we use a sixth-order 
polynomial p to fit the correction term. The figure shows that the fit is significantly worse than 
in figure Rl See the text for a discussion of possible reasons for this effect. 
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shock (see e.g. [CIR52, MRTB05]), but even in this case we expect the shock to propagate at 
the wrong speed. 

Another problem is the following one: In the case of a non-conservative discretisation of 
the type f|33|) . a simple linear stability analysis reveals that the discretised solutions develop 
an ultraviolet instability (i.e. the mode v n — (—1)™ becomes unstable) in the regions where 
u > 0. Similarly, the corresponding left-sided discretisation shows an instability in the regions 
with u < 0. In contrast, for the case of a centred discretisation, the highly oscillatory modes are 
stable independently of the sign of u. 



Remark 4.1 Conservative discretisations, i. 

1 
2d 



d t U n = —^(ul +1 -ul) 



or variants thereof, still require the scheme to be upwind for stability, but in this case the shocks 
of the discretised system propagate at the correct speed. Also, these schemes do not suffer from 
the ultraviolet instability. 

How is this picture modified for non-zero values of el The instabilities discussed above grow 
at a speed 0(S^ 1 ) while the stabilising effect of the viscosity is of the order e5~ 2 . Thus, one 
expects the viscosity to dominate only if e ^> 6. Regarding the behaviour after the formation of 
a shock, a simple boundary layer analysis shows that a typical shock for ([52)1 has width O(e), 
so that the caveats pointed out above are expected to become relevant as soon as e < 8 (see for 
example EVEOOa for a more sophisticated boundary layer analysis that even goes to the next 
order in e). 

On the other hand, at least away from shocks, the analysis performed in section [2] holds as 
soon as u can be approximated by the solution to the linearised equation at sufficiently small (but 
still much larger than 5) spatial scales. The fcth mode presents spatial features on a lenghtscale 
of order fc _1 and the nonlinearity essentially propagates the solution through space at speeds 
of order 1. Thus, the timescale on which the fcth mode varies due to the nonlinear effects is 
expected to be of order fc _1 . On the other hand, the relevant timescale of the linear part for 
this mode is (efc 2 ) -1 and we can conclude from this heuristic consideration that the linearised 
equation is a good approximation to the full solution for modes with fc ^> s , i.e. on space 
scales much smaller than e. Consequently, one again expects the results from section [5] to be 
relevant as long as e ^> S. This leads to the following statements. 

Conjecture 6 For e«1, we expect the solution to the finite difference approximation of A Sty) 
to show the following behaviour. 

1. For 8 <^i e, the discretised solution converges to the viscosity solution of 

d t u = ~d x (u 2 ) + ~, (34) 

where c € {1,0,-1} depending on whether the discretisation is right-sided, centred, or 
left-sided. Neither of the instabilities discussed above occur. 

2. For £<5, both viscosity and the noise term become irrelevant; the solution behaves like 
the corresponding approximation to the inviscid Burgers equation. In particular, as long as 
there are no shocks and while solutions have the correct sign to prevent ultraviolet blow-up, 
we expect to converge to \3m with c — 0. After the occurrence of a shock, one expects 
stability only if the scheme is upwind and, for discretisations of type 133\) . shocks will have 
the wrong propagation speed. 

To test this conjecture we again perform a numerical experiment. Keeping in line with the 
topic of this article, we focus on studying how the presence of the extra c/4 term in (|3"4"]) is affected 
by the values e, 5 > 0. For our experiment, we first solve the limiting equation (1341) numerically 
up to a time t > 0, once with c = and once with c = 1, to get states Uo(t), «i(£) G L 2 ([0, 27r],R). 
Now, for given e, 5 > 0, we solve the right-sided finite difference discretisation (fT3"|) for ([3"2"| to 
get a solution u e< s(t). According to conjecture El we expect u e> s{t) ~ U\ for 5 <C e <C 1 and 
u e ,s{t) ~ u o f° r £ "C S. 
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Figure 7: Illustration of the limiting behaviour of a right-sided finite difference discretisation 
for 1 32(1 as e \. (7or /ixeci 5J. Distances are shown in "two-centre bipolar coordinates" as 
described in the text. One can see that, as e gets small, the finite difference discretisation first 
gets close to the solution of \3J$ for c — 1 but then finally converges to the solution with c = as e 
gets smaller than S. This behaviour reflects the dichotomy between the two cases of conjecture® 



To verify this conjecture, in figure[7]we consider, for fixed S, the solution u Sj s{t) as a function 
of e. The coordinate system is chosen such that the (two dimensional) distance of u e ,s(t) to 
the point c = in the graph equals kt £]( 5(i) — ito|| i 2; similarly the distance of u et s{t) to the 
point c = 1 in the graph equals | kt E)1 5(t) — Mi||i,2 (this is sometimes called "two-centre bipolar 
coordinates"). The simulation indicates that the transition from c = 1 to c — takes place at 
around e « S as expected. 

A Simulations 



To verify the heuristic arguments presented above, the text utilises a series of numerical results. 
This appendix summarises some of the technical aspects of these simulations! 3 ] 

We first describe how to implement the finite difference schemes used to discretise SPDEs 
like (HJ, (J2TJ) and (|29)): For the space discretisation of these equations we approximate states 



G L 2 ([0,2tt],R) by vectors u 



N 



pJV 



The space discretisation of the differential operators 



given by formula (|T^|) . above. The finite difference discretisation of the white noise process w is 
W N /y/Ax, where W N is a standard Brownian motion on R N and Ax = S = 2n/N is the space 
grid size. This leads to R^-valued stochastic differential equations of the form 



du N = vL N u N dt + F N (u N )dt + o(u N ) \ -^dW N (t) 

" Ax 



where Ln £ M. x is the discretisation of the linear part and Fn : 
of the nonlinearity. 

For discretising time we use the 0-method 



dN 



pN 



is the discretisation 



u (n+i) = u (n) + i ,L JV (6» u («+i) + (i _ 0)u M ) At 

+ F N (u^)At + a(u^)^e n \ 

V Aa; 

where At > is the time step size, u^ is the discretised solution at time n At, the ^ n - ) are i.i.d., 
A^-dimensional standard normally distributed random variables, and 6 = 1/2. Rearranging this 



3 For full details we refer to the source code of the programs used in these simulations, which is available for 
download at http://seehuhn.de/programs/HairerVosslO . 



16 



equation gives 

(I - v6AtL N )u^ n+1) = (I+ i/(l - 6)AtL N )u^ 

n^t (35) 

+ Fjv(w (n) ) At + a(«(»))«/_lf(»). 

V Aa; 

Relation (l3"5j) allows to compute u tjl+1 ^ > from u' n ); since / — vOAtL^ is cyclic tridiagonal, this 
system can be solved efficiently. 

Since we are using the partially implicit 0-method for the linear part of the SDE, there are 
no constraints on the time step size At arising from this term; on the other hand, since the 
non-linear transport term is treated explicitly, one has to make sure that vAt/Ax < C, where 
v is the largest speed of propagation appearing in the solution and C is the Courant number, 
thus ensuring that the CFL condition is satisfied. The resulting method can be used to perform 
the simulations required to generate figures [TJ [5] and [3 as well as the "finite difference" curve in 
figured 

As described in section [2J we can compute a spectral Galerkin approximation to the second 
order differential operator using discrete Fourier transform. This corresponds to replacing Ln 
in (|35[) with 

L N = T~ X DT 

where D = diag(-0 2 , -l 2 , . . . , -|_iV/2j 2 ) and J 7 : R N -> C^ N/2i+1 represents the discrete Fourier 
transform (since the data is real, only L^/2J +1 of the Fourier coefficients need to be considered). 
Because the matrix Ljv is no longer tridiagonal, one should not try to explicitly construct this 
matrix. Instead one can use the fact that T and T~ x can be computed efficiently: we can 
compute the right-hand side of ((35]) using 



(I + i/(l - 6)AtL N )u (n) = T~ x diagfl - v{\ - 6) At k 2 , k = 0,..., [N/2\ 2 ) Fu (n) 
and to solve (|55jl for u^ n+1 ^ we can use the relation 

(I veAtL N )-h = J- 1 diag( i + ^ Atfc2 , k = 0, . . . , [N/2\ 2 )Tb. 

This technique allows to obtain the "Galerkin" curve in figure [2J Similarly, the rougher-than- 
white noise for figure [3] was obtained by replacing the noise term £( n ' with 

£(«) = J- 1 diag((l + fc 2 )-' 5 , fc = 0,...,LA r /2j)J"C ( ™ ) - 

Finally, for the minimisation procedure performed to generate figures [4] and [6] we employ the 
simplex algorithm by Nelder and Mead NM65 . 
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